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Abstract. - The behaviour of the one-dimensional random-forced Burgers equation is investi- 
gated in the path integral formalism, using a discrete space-time lattice. We show that by means 
of Monte Carlo methods one may evaluate observables, such as structure functions, as ensemble 
averages over different field realizations. The regularization of shock solutions to the zero-viscosity 
limit (Hopf-eq.) eventually leads to constraints on lattice parameters, required for the stability of 
the simulations. Insight into the formation of localized structures (shocks) and their dynamics is 
obtained. 



I. Introduction. — The random-force-driven Burg- 
ers equation 

dtu + udxU — vd^u = f (1) 

itiay be interpreted as a model for compressible hydrody- 
namic turbulence, describing acoustic perturbations in the 
I'eference frame moving with the velocity of sound [1,2]. 
Here u is the velocity, / the random forcing, and v is the 
kinematic viscosity. One generally assumes the force to be 
white-in-timc, displaying Gaussian statistics 



(2) 



so that the properties of the external pumping arc com- 
pletely characterized by the covariance x- The turbulent 
state is maintained by a large-scale external force with 
correlation length L. Finite i> > and energy dissipa- 
tion e = x(0) provide a dissipation scale ?7 = (i^^/e)^/'' > 0. 
The dissipation is related to the Reynolds number Re by 
the definition e = Kc^v^/L^. Thus the characteristic ve- 
locity is ul ~ (eL)^/''. In the limit of large Reynolds 
number Rc = (L/rj)^/^ ^ 1 these scales separate, so one 
expects a turbulent cascade. Evaluating moments of 
Galilean invariant velocity differences w{r) = u{r) — u{0), 
i.e. structure functions 



5,(r) ^ {[w{rW), 



(3) 



p > 0, where r is the displacement, one therefore expects a 
scaling behavior in the inertial interval L ^ r ^ rj. From 
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dimensional considerations we find Sp{r) oc (er)'>p, and 
scale-invariance implies Cp ~ p/S- However, the formation 
of shocks in Burgers turbulence leads to a strong intermit- 
tency and a bifractal scaling of the form (^p = min(p, 1), 
see e.g. [3]. 

We would like to emphasize that the average ( • • • ) in 
Eq. (3) and in the following is usually considered to be an 
average over all times at a given place, or over all places 
at a given time; linked by the Taylor-hypothesis in sta- 
tistically homogenous turbulence. In this article these av- 
erages are represented as ensemble averages, as discussed 
below. 

Numerical investigations of hydrodynamic systems are 
commonly performed by means of direct simulations or 
variants thereof which amounts to integrating the equa- 
tions of motion. In this article we follow a completely 
different approach, namely to evaluate the correspond- 
ing path integrals by Monte Carlo simulations. The clear 
physical picture of Burgers equation, and the fact that the 
intermittent structures are well-known - they correspond 
to shocks with a large negative velocity gradient - makes 
the Burgers equation an attractive setting to study Monte 
Carlo methods in turbulence theory. Further, a number 
of technical advantages leads to the unique possibility of 
investigating intermittent structures in turbulent flow. In 
particular. Burgers equation is local while incompressibil- 
ity acts as a nonlocal interaction in Navier-Stokes turbu- 
lence which eventually complicates issues due the presence 
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of the functional determinant arising in the path integral 
formulation of Eq. (1) (see Sec. II). Also the dissipation 
scale rj provides a UV-regularization of shock structures. 
Finally, a huge variety of analytical methods has been ap- 
plied to Burgers equation, giving results that can directly 
be compared to numerical calculations. However, Burgers 
equation is certainly also interesting on its own (e.g. in 
cosmology), serving as model for systems with an inter- 
play of random processes and coherent structures; for an 
overview see [3] . 

It is the aim of this paper to demonstrate the feasibility 
of Monte Carlo simulations in the path integral formal- 
ism. In particular we will show that respecting certain 
constraints on the parameters of the model, stable simu- 
lations can be performed. We will demonstrate with an 
example that this way structure functions can be com- 
puted precisely. 

This paper proceeds as follows. In Sec. II we shortly 
introduce the path integral formulation for Burgers equa- 
tion. Then in Sec. Ill the numerical scheme is discussed 
in detail and in Sec. IV the numerical results are given. 

II. Path Integral for Burgers equation. The 

path integral for the randomly forced Burgers equation 
is introduced via the Martin-Siggia-Rose formalism [4] 
where one has 



Then the expectation values of Galilean invariant ob- 
servables 0[m(x)] are calculated as ensemble-averages 



Z= DuDp cxp{- S[u,p]) 



(4) 



with the conjugate momenta p, and the action is 
S[u,p] ~ Jdt C[u,p]. The Lagrangian is given by 



C = —i J dx p (^dtu + udxU — vd^u) + — Jdxp{x*p), (5) 

where * denotes the convolution. Performing the Gaussian 
integration over the conjugate field, one arrives at 



(6) 



where N[u] = dtu + udxU — v&^u. The functional deter- 
minant, arising in principle in the derivation of (4), (5) 
does not contribute to the forced Burgers equation due to 
the local interaction (see e.g. [5]). Thus the integration 
measure is properly defined, and the integration is over all 
fields u (and p) in one dimension. In contrast, in the case 
of the Navier-Stokes equation, where the effect of pressure 
is not negligible, the corresponding non-local interaction 
leads to a non-vanishing functional determinant that has 
to be treated using the Fadeev-Popov procedure [6] . 

As (4), (5) are Galilean invariant, the computation of 
non-invariant quantities such as n-point correlation func- 
tions would require a gauge fixing [7]. However, since we 
are interested in evaluating structure functions, i.e. man- 
ifestly Galilean-invariant quantities, this procedure does 
not have to be performed. 



(ON) = 



jDuO[u] exp(-S'N) 
/Du cxp(-S'N) 



(7) 



Extending the time integration in the action to infinity 
renders the averages stationary (see e.g. [8]). 

Using the path integral (4), (5) it has been shown that 
intermittent statistics of Burgers equation can be under- 
stood in terms of instanton solutions [9-11]. These solu- 
tions break Galilean invariance which has been shown to 
account for intermittency [12]. 

III. Monte Carlo Simulations. Using the path 
integral formalism we calculate ensemble averages like (7) 
numerically, where the Lagrangian is given by (6) . Defin- 
ing the theory on a finite size lattice, the integration is 
performed over grid variables Ui(n), i.e. the field u defined 
at lattice sites (i, n) where i denotes the discretized space, 
and n discretized time points. Then we obtain 



(O) = j\[du,{n) 0[{u,{n)}] exp ( - 5[K(r^)}] 



(8) 



for the observable computed in the discretized theory. 
Here the discretized action S[{ui{ri)}] may be so chosen 
that the single-site action S(ui{n)^ is quadratic: 



S{u,{n)) = A,{n)[u^{n) - C^{n)Y + B,{n) 



(9) 



where Ai{Ti) > 0. The terms Ai{n), Ci{n) and Bi(n) 
will in general depend on the field at different nodes on 
the lattice. The action is locally minimized by choosing 
Ui{n) ~ Ci{n) thus giving the main contribution to the 
integral in (8). This way one may evaluate the whole inte- 
gral, by successively integrating over different lattice sites. 
Successive application of such sweeps, will then generate 
the Markov chain {uf^\n)} of configurations in the fc-th 
step of the process. The quadratic form of the action in 
(9) allows to apply a heat-bath Monte Carlo algorithm to 
update single field variables on the lattice. However, it 
should be emphasized that writing the single-site action 
as (9) imposes a restriction on the possible form of the 
nonlinear term udxU. In its discretized version, it has to 
be at most linear in Ui{n). Thus we are led to define the 
nonlinear term as udxU — *■ Ui{n)u[{n) , where the deriva- 
tive is written symmetrically with the simplest choice be- 
ing u'^{n) = (ui+i{n) — 7ii_i(n))/2Ax. The time deriva- 
tive and the viscous term are discretized in the usual way. 

The above action was put onto a rectangular lattice with 
A'x sites in space- and A^^; sites in time-direction. Lat- 
tice spacings are denoted correspondingly by Ax and At. 
Monte Carlo simulations were performed for different lat- 
tice sizes, ranging from e.g. (A'x = 16) x (A^^^ = 4), to larger 
lattices of up to approximately (Ax = 362) x [N^ = 4096) 
at fixed viscosity respectively. The Reynolds number 
ranged from Re = 2 to Re = 16. 
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Lattice parameters and Continuum Limit. In our sim- 
ulations we chose the forcing characteristics as 

xix-x') = xiO)ei^p{-\x- x'\/L) , (10) 

in agreement with the considerations in [10]. In Fourier- 
space this represents a noise x(fc) with an IR-cutoff at 
wave- number k ~ L~^ and a power-law behavior for 
k > L-i, of the form x(0) A;^, where /3 = -2. That 
is a random force acting predominantly at large scales > L, 
leading to shock-formation at random positions in our sys- 
tem [3]. 

To identify the lattice parameters with the correspond- 
ing constants of the continuum theory, first notice that the 
viscosity has to be defined as 

1^ = i){Ax)^/At (11) 

where the arbitrary constant i> has been introduced, relat- 
ing units for space- and time-measurements on the lattice. 
This can be shown by performing the continuum limit in 
the case of the symmetric random walk, leading to the dif- 
fusion equation. The so-defined ly gives us the Reynolds- 
number Re according to the earlier given relations, after 
x(0) and L have been specified. The continuum limit is 
then performed by Ax, At — > while keeping L, v and Re 
constant. 

In both space- and time-direction periodic boundary 
conditions were specified. Then one may expect that in 
the continuum limit the correct physical behavior of the 
fields is assumed and correlation effects induced by the 
boundary conditions are negligible. 

Stability considerations lead to further constraints for 
the lattice size, as will be explained in the following sub- 
section. 

Stability and Regularization. As one would expect, the 
stability of the simulations over a large number of Monte 
Carlo steps is a big issue. Indeed, if certain restrictions on 
the lattice parameters are not taken care of, the simulation 
terminates sooner or later due to instabilities. 

The occurence of instabilities in our Monte Carlo sim- 
ulations is related to the existence of shock-like solutions 
of Eq. (1). To obtain stable simulations, the shock struc- 
tures have to be resolved on the lattice, otherwise the ac- 
tion near the shock positions is represented badly and the 
simulation tends to develop instabilities, accumulating the 
overall energy of the configurations near the shocks. 

In the limit of small viscosities it would be possible to 
regularize the shocks by adding higher derivative terms 
or nonlinear terms to the differential equation. Simi- 
lar to the case of the Korteweg-de Vries equation, where 
the term proportional to d^u serves to regularize shock 
waves [13, 14], the Burgers equation can be regularized 
by adding e.g. terms proportional to dtd^u and ud^u to 
it [15]. This would smear out the shocks and enable their 
proper treatment in the simulation. As indicated by the 
shock solution of the unforced Burgers equation 

u = -Atanh —x , (12) 



the viscosity term itself provides a valid regularization of 
shocks. Therefore we decided to smooth out shocks by 
chosing the viscosity appropriately. This amounts to re- 
solve the dissipation scale 77 on the lattice. In terms of the 
Reynolds number this translates into the constraint 

Re<LiVx. (13) 

At given Reynolds number Re 3> 1 and L < 1 this imposes 
a strong restriction on the possible lattice sizes. 

Respecting the constraints on the parameters we are 
in the position to perform stable simulations. In partic- 
ular, we can resolve a shock on our lattices as Fig. (1) 
demonstrates, which shows a time-slice of a typical con- 
figuration. We also can observe time dynamics in our sim- 
ulations as Fig. (2) shows, and finally, we are able to com- 
pute the structure functions an example of which is given 
in Fig. (3). 

IV. Summary and Outlook. — We have shown how 
to perform stable Monte Carlo simulations of stochastic 
partial differential equations like the Burgers equation in 
the path integral formulation, and in this way demon- 
strated that it is possible to calculate observables like the 
structure functions using an ensemble-average rather than 
a time- or space-average. The lattice versions of the the- 
ories can directly be identified with their continuum coun- 
terparts, and as long as certain constraints on the lattice 
size are respected no instabilities occur. Direct insight into 
the structures leading to intermittency and, thus, multi- 
scaling, can be obtained. Especially, we want to point 
out that the existence of a dissipation length scale can be 
observed. 

In particular in our simulations of the randomly forced 
Burgers equation we observe the following: 

• Thermalization and autocorrelation times are very 
long, up to the order of 10^ Monte Carlo steps. 

• In the stable runs after thermalization, the typical 
shock solutions of Burgers equation and their dynam- 
ics can be observed, see Fig. (1) and Fig. (2). 

• In the unstable runs before occurence of instabilities, 
configurations show typical shock solutions leading 
eventually to an energy accumulation near shock po- 
sitions. 

• The distinction between stable and unstable simula- 
tions can directly be related to the existence of a dis- 
sipation length scale which is either bigger (stable) or 
smaller (unstable simulations) than the lattice spac- 
ing. 

First results on the extracted structure functions, e.g. see 
Fig. (3) are promising. However, the understanding of 
systematic errors such as finite size effects or fitting ranges 
for extracting the scaling exponents from the measured 
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Fig. 1: Slize of a typical configuration u{x,to) at constant time 
to, depending on the (spatial) lattice site x/Ax; taken from a 
iVx = 256, = 45 lattice where Re = 4, i/ = 1/32. The typ- 
ical kink solution associated with the Hopf-eq. can be clearly 
seen, including a finite shock width due to regularizing viscos- 
ity. 
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Fig. 2; 3D-Plot of a configuration u{x,t); taken from a TVx = 
32, A^^^ — 64 lattice. The propagation of the shock with con- 
stant velocity is clearly visible. 



structure functions is still in progress. Nevertheless, we 
conclude that the path integral Monte Carlo approach to 
the Burgers equation outlined in this article is a feasible 
and stable tool to obtain interesting numerical results. 



Fig. 3: Third order structure function 53 (r) as a function of 
the separation r/Ax, evaluated from a simulation with 5 * lO'' 
configurations on a A'x ~ 256, A'^ — 45 lattice with Re = 
16,1^ = 1/512. 
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